echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 9.S (Runge-Kutta Method for a higher order).
% Section	9.7,	Runge-Kutta Methods, Page 475
echo on; clc; format long; hold off; clear

% This program implements the Runge-Kutta method
% for solving a second-order initial value problem:
%     x'' = f(t,x,x')     with     x(a) = xa  
%                                 x'(a) = ya

% Which is equivalent to solving the system:
%      x' = y             with     x(a) = xa
%      y' = f(t,x,y)               y(a) = ya

% The formula for f(t,x,y) and is used to form the
%	vector function F(t,Y)	which is stored in  fn.m 

% function W = fn(t,Z)
% x = Z(1);  y = Z(2);
% W = [y, -40*sin(x)];

delete fn.m
diary fn.m; disp('function W = fn(t,Z)');...
            disp('x = Z(1);  y = Z(2);');...
            disp('W = [y, -40*sin(x)];');...
diary off;

fn(0,[0 0]); % Remark.  fn.m and rks4 are used in Algorithm 9.H
pause	       % Press any key to continue.

clc;
%	Place the endpoints of [a,b] in  a  and  b.

%	Place the initial value  Z(a) = [xa ya] in  Za.

%	Place the number of subintervals in  m.

a  = 0;
b  = 2;
Za = [0.3  0.0];
m  = 40;
[T,Z] = rks4('fn',a,b,Za,m);
P = [T;Z(:,1)']';
points = P(1:4:length(P),:);

pause	% Press any key to see the list of solution points.

clc;, clg;
X = Z(:,1);
a = min(T)-0.3;
b = max(T)+0.3;
c = min(X)-0.3;
d = max(X)+0.3;
axis([0 2 -0.35 0.35]);...
plot(T,X,'g');...
hold on;...
plot([0 2],[0 0],'r',[0 0],[-0.35 0.35],'r');...
xlabel('t');...
ylabel('x');...
title('Runge-Kutta solution to x`` = f(t,x,x`)');...
grid;...
axis;...
hold off;...
shg; pause	% Press any key to continue.

Mx1 = 'Runge-Kutta solution to x`` = f(t,x,x`).';
Mx2 = '     t(k)               x(k)';
clc,echo off,diary output,...
disp(''),disp(Mx1),...
disp(''),disp(Mx2),disp(points),diary off,echo on
